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ABSTRACT 

We present a detailed study of the hydrodynamical interactions in a 
Newtonian black hole-neutron star binary during the last stages of inspiral. 

We consider close binaries which are tidally locked, use a stiff equation of state 
(with an adiabatic index P = 3) throughout, and explore the effect of different 
initial mass ratios on the evolution of the system. We calculate the gravitational 
radiation signal in the quadrupole approximation. Our calculations are carried 
out using a Smooth Particle Hydrodynamics (SPH) code. 

Subject headings: binaries: close — gamma rays: bursts — hydrodynamics — stars: 
neutron 


1. Introduction 

It is our objective in the present work to investigate the outcome of the coalescence 
of a black hole with a neutron star, to hnd out to what extent the neutron star is tidally 
disrupted, and in particular to determine if an accretion structure does form around the 
black hole as a result of the encounter. We also wish to explore the role of the mass ratio 
and the stiffness of the equation of state in the evolution of the system and in the emission 
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of the gravitational radiation. For simplicity and to perform accurate comparisons with the 
work presented by [Rasio fc Shapiro 1994| , (hereafter RS) for double neutron star binaries, 
we have chosen to model the neutron star as a polytrope and explore different mass-radius 
relationships by varying the stiffness through the adiabatic exponent F. In this paper we 
restrict ourselves to a stiff equation of state with F=3. 

We restrict our work to binaries which are tidally locked. The construction of a 
completely self-consistent initial condition is then entirely straightforward. This assumption 
represents a major simplihcation in our calculations, and can be seen as an extreme 
case of angular momentum distribution in the system. The opposite extreme would be a 
case where the binary components exhibit no rotation at all as viewed from an external, 
inertial frame of reference, and it poses a much greater problem regarding the construction 
of self-consistent initial conditions in equilibrium. Since complete tidal locking is not 
expected ( [Bildsten &: Cutler 1992|) , the full range of initial conhgurations deserves to be 
explored. 

The motivation for our work is presented in section The numerical method we have 
used to carry out our simulations is described in section Calibration to previous work and 
our modeling of the black hole-neutron star system is presented in section In section ^ 
we present our numerical results, and a discussion of these results follows in section p. 


2. Motivation 


The orbital evolution of short-period binary systems containing compact objects (black 
holes or neutron stars) will inevitably lead to coalescence due to orbital angular momentum 
loss to gravitational radiation. Such binary systems have been found (PSR 1913+16— [Hulse 
I&: Taylor 1975| , PSR 1534+12— [Wolszczan 199 1|) and the rate of inspiral matches the 
prediction of general relativity to high accuracy ([Taylor fc Weisberg 1982| , [Taylor et al. 
|1992|) . The time to hnal merging in these binaries is less than the Hubble time and so 
they will eventually coalesce. In the hnal stages of such a coalescence, a powerful burst of 
gravitational waves is expected, and as such these systems are primary candidate sources 
for detection by instruments such as LIGO and VIRGO, expected to begin operation 
within a few years ( [Fhorne 1995| ). Estimates for the event rate can be inferred from the 
statistics of the known Hulse-Taylor type binaries, and a rate of 10“® to 10“® per galaxy per 
year can be expected ([Lattimer fc Schramm 1976[ [Narayan, Piran fc Shemi 1991|) , which 
would imply several coalescences per year could be observed out to a distance of 1 Gpc. 
These gravitational waves will certainly carry a copious amount of information about the 
emitting source. Among other things, they are expected to exhibit strong departure from 
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point-mass behavior due to hydrodynamical effects and the ffnite size of the stars, as well 
as a dependence on the details of the equation of state for matter at nuclear densities. 

A fully three-dimensional hydrodynamical study is therefore essential to understand 
the gravitational radiation signal, and this information can be used to place constraints 
on the equation of state. For the typical masses and binary separations involved in very 
close binary systems (a few solar masses and several tens of kilometers respectively), one 
should ideally solve the problem using general-relativistic formalism and include radiation 
processes. This has not yet been realized, but some qualitative results can hopefully be 
obtained from a Newtonian treatment. 


It has been shown (|Lai, Rasio fc Shapiro 1993a|) that at least in the Newtonian 
case, hydrodynamic effects play an important role in the orbital evolution of the system. 
Essentially, tidal interactions can make a close binary dynamically unstable for small enough 
separations—on the order of a few stellar radii—when the spin angular momentum becomes 
comparable to the orbital angular momentum. This effect alone can produce orbital decay 
on a time scale comparable to that of angular momentum loss to gravitational waves. 

The final stages of inspiral, and the coalescence of double neutron star binaries 
has been the object of numerous previous studies in a Newtonian framework. The 
gravitational radiation signal was studied by Nakamura Sz Oohara (|1989|) and Oohara & 
Nakamura (|1989| , |1990|) while Davies et al. (|1994|) considered the effect of different initial 
spin configurations on the outcome of the merger. More recent simulations have calculated 
the neutrino emission as well as the gravitational waves from the coalescence (Ruffert et 
al. |1996| , |1997a| , |1997b|) . The de-stabilization effect of tidal forces in a Newtonian framework 
has been investigated in great detail (|Lai, Rasio fc Shapiro 199?5| , RS) using an analytical 
approach based on an energy variational method to examine equilibrium configurations 
and a three-dimensional numerical treatment to study the hydrodynamical aspects of the 
coalescence, using throughout a polytropic equation of state. 

There is yet another reason why these systems are the focus of intense study, since they 
have been suggested ( [Eichler, Livio, Piran fc Schramm 1988|) as a possible source for the 
production of gamma-ray bursts (GRBs). It is well established that the GRB distribution 
on the sky is isotropic and it is generally believed that their sources lie at cosmological 
distances ([Fishman fc Meegan 19951) , with redshifts on the order of unity ( [Metzger et al.\ 
1997[) . The observed ffuence of said bursts, 10“^'^ to 10“^ erg cm“^, and the rise times 
in their light curves imply that an energy release of at least 10^^ erg takes place in a few 
milliseconds, in a region at most a few hundred kilometers across. 


To produce a GRB, this amount of energy, after being released in a primary event. 
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must be transformed into 7 -rays. The relativistic blast wave model ([Meszaros fc Re^ 

requires a relatively baryon-free line of sight to the observer along which matter 
can be accelerated to speeds close to the speed of light. The interaction of this outflow with 
the interstellar medium would then presumably lead to shock acceleration of electrons and 
emission of 7 -rays through synchrotron radiation. 


1992,1993 


In the double neutron star merger scenario (Paczyhski 1986| , [Eichler, Livio, Piran 
fc Schramm 1988|) the merger initially produces a burst of neutrinos and anti-neutrinos. 
However, Newtonian calculations (|Janka &: Ruffert 1996|) have shown that the energy 
release into neutrinos is insufficient to power the observed GRBs. General relativistic 
calculations carried out in a conformally flat metric (|Wilson, Mathews fc Maronetti 1996| ) 
have shown that in some cases each neutron star may collapse into a black hole many orbits 
prior to merging so that no blast wave will occur. Nevertheless, this could still produce 
a GRB with a smooth time prohle ([Wilson, Salmonson fc Mathews 1997|) . In the black 
hole-neutron star merger scenario (Paczyhski 1991|) , it was expected that the neutron star 
will be disrupted by the black hole and a thick accretion torus will be formed. As in the 
case with two neutron stars, this could lead to the formation of a blast wave which would 
power the GRB. 


3. Numerical Method 


For the calculations presented in this paper, we have used the numerical technique 
known as Smooth Particle Hydrodynamics (SPH). This is essentially a Lagrangian method, 
where forces are evaluated by interpolation over a grid of points co-moving with the fluid, 
which can be considered as particles. This method was originally developed by Pucyj 
1(1977)1 and [Gingold fc Monaghan (1977)| as an alternative to Eulerian hydrodynamics 
(computations on a hxed grid). The principal advantages of SPH are that no assumptions 
need to be made a priori about the nature of the flow that will be studied and that no 
computational effort is wasted in modeling regions where matter is not present. This 
is particularly advantageous when studying complicated three-dimensional astrophysical 
flows. An excellent review of the method has been given by [Monaghan (1992']| . 


SPH has been tested successfully and applied to a variety of problems, such as the 
shock wave in a tube problem ( [Monaghan fc Gingold 1955[ ), redistribution of angular 
momentum in a thick accretion torus ([Zurek fc Benz 1986[) , static stellar structure ( pingold 
&; Monaghan (1977][) , astrophysical jets ([Goleman fc Bicknell 1985[) , and hydrodynamics of 
close binary systems ( Penz et al. 1990 [ , [Rasio fc Shapiro 1992[ , [1994[ ). We have developed 
our own SPH code ([Lee 1998[) and successfully tested it in one, two and three dimensions 
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with several of these problems. Agreement with previous results has been excellent in every 
case. A description of our code is given in Appendix 


4. Calibration and initial conditions 


4.1. Simulation of a double neutron star binary 


Given that the problem of merging black hole-neutron star binaries is closely related 
to the study of coalescing double neutron star or white dwarf binaries, we have performed 
a rigourous calibration of our code to the results presented by RS for two polytropes with 
stiff equations of state (F=3) and an initial mass ratio of unity in a tidally locked binary. 
While our simulations had a lower resolution (~2000 particles per star vs. 40000 for RS), 
qualitative and quantitative agreement was excellent. 

For the following, measure distance and mass in units of the radius and mass of the 
unperturbed (spherical) neutron star (13.4 km and 1.4M0 respectively), except where 
noted, so that the units of time, density and velocity are: 


i = 1.146 X 10"^ s X 


R 

13.4 km 


3/2 / ^ \ 1/2 

[ti^) 


( 1 ) 


P 


1.14 X 10^® kg m ® X 


R 


13.4 km 


-3 


M \ 
IAMq ) 


( 2 ) 


V = 0.39c X 


R 


13.4 km 


- 1/2 


M 


1/2 


1.4Mq 


( 3 ) 


where R and M are the radius and mass of the unperturbed (spherical) neutron star. 


The coalescence resulted in a massive central core containing about 82% of the total 
mass, surrounded by a massive halo and extended spiral arms. At the end of the calculation, 
the central object is not azimuthally symmetric and thus continues to emit gravitational 
waves, of amplitude hfi^ai- The amplitudes and frequencies of the gravitational radiation 


waveforms were found to be in excellent agreement with RS (see |Lee fc Kluzniak 1995|, [Lee 
19981) . In Table |l| we present a comparison of some of the more important parameters in the 
case of the coalescence of two neutron stars. The first, second and third columns show the 
maximum and hnal amplitude in the gravitational radiation waveforms and the maximum 
gravitational wave luminosity respectively. These values are given in geometrized units. 












6 


where G=c=l, and the peak luminosity Lmax is normalized to Lo=c^/G=3.59 x 10^® erg s“^. 
The fourth and fifth columns display the central density of the resulting core and the mass 
contained in the halo surrounding it, as a fraction of the total mass in the system. 

For the simulations of a black hole-neutron star binary presented here, we have chosen 
to keep a stiff equation of state (r=3) to compare our results to the case of the coalescence 
of two neutron stars. In future work, we will explore different values of the adiabatic index 
r to model different mass-radius relationships. This is particularly important since future 
measurements of the gravitational radiation emitted during these events could serve to 
constrain the equation of state for dense matter. 


4.2. Modeling of the black hole 


In our simulations, the black hole is modeled by a point mass, and produces a 
Newtonian potential: 




GMbh 

r 


The contribution from the black hole, to the force on particle i in the star is then simply 
given by: 


pBH 

i 


GMsH^i 

ri - vbh 




f'BH), 


and symmetrically. 


^ BH ~ 


GrriiM, 


BH 


Ti — Vbh 


[rBH 


n), 


is the contribution from particle i to the force on the black hole. 

The horizon of the black hole is modeled by placing an absorbing boundary at a 
distance rs=2GMBH /from the point mass. At every time step during the dynamical 
simulations, any particle that crosses this boundary is absorbed by the black hole and 
removed from the simulation. The mass, position and velocity of the black hole are adjusted 
so that total mass and linear momentum are conserved. We thus disregard any spin angular 
momentum that might be gained by the black hole during the process. This does not 
present any problems, since our calculation is Newtonian throughout, and we make no 
attempt to model frame-dragging. 
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4.3. Construction of initial conditions 

Our simulations begin with the construction of a spherical, unperturbed polytrope, 
which for simplicity we will refer to as the neutron star. To do this, we proceed essentially 
as Rasio & Shapiro (1992). A total of N particles are placed on a cubic lattice with masses 
mi=pLE/n, where ple is the density calculated from the Lane-Emden solution to the 
equation of hydrostatic equilibrium with an equation of state P=Kp^ and n is the number 
density of particles. The smoothing length is assigned to each particle so that the number of 
overlapping neighbors per particle is z/=64. This amounts to setting the smoothing length 
hi ~ /, where I is the lattice spacing. There is considerable advantage in having variable 
particle masses, as this increases the spatial and density resolution near the edge of the star, 
where the density gradient is largest. The radius of the unperturbed polytrope is 13.4 km, 
and its mass is 1.4M0. For the simulations presented in this paper, we have used N ~ 8000 
for every case except two, where N ~ 17000. 

To obtain tidally locked equilibrium conhgurations for the binary system, the 
unperturbed polytrope and the black hole are placed in the co-rotating Keplerian frame of 
initial angular velocity f2= ijG{M + Mbh)/^^, where r is the binary separation dehned as 
the distance between the black hole and the center of mass of the neutron star. A damping 
term linear in the velocity is introduced into the equations of motion for the SPH particles 
to allow the star to respond to the presence of the tidal gravitational held. We now have: 

Fiq T Fifj '^i/tdamp T Efli^l 

Mbh'^bh = FbhPMbh^'^'t'bh- 

where FiG and FiH are the gravitational and hydrodynamical forces on particle i (of mass 
rrii) respectively. The value of tdamp is chosen so that oscillations are critically damped, 
reaching an equilibrium as quickly as possible; tdamp = {GM/R^)~^. We neglect Coriolis 
forces since we are interested in equilibrium conhgurations with no bulk motion in the 
co-rotating frame. The position of the black hole and the center of mass of the star are 
adjusted at every time step so that the binary separation remains at a desired value. In the 
same manner, is adjusted so that the total (gravitational plus centrifugal) force on the 
black hole in the co-rotating frame, FEH+MBH^^rSH is zero, i.e.: 


where vbh is the position vector of the black hole. This procedure ensures that the 
conhguration reaches equilibrium in a state of synchronization. 




For a given value of the binary separation, we allow the system to relax for a period of 
twenty time units, keeping the specific entropies of all particles constant, i.e. iF=constant 
in P=Kp^. In all cases, our initial conditions satisfy the virial ratio to better than three 
parts in 10^. We have varied the mass ratio in the binary, q=M/M bh, by changing the 
mass of the black hole only. 


5. Results 
5.1. Introduction 

In this section we present the results of our simulations for different values of the mass 
ratio q in the binary, beginning with high values of q and proceeding in descending order. 
The units used are as dehned in equations (|l|), (|D and (|^) except where noted. The highest 
value of the mass ratio is unity. Although the production of a binary system with such a 
low-mass black hole is rather unlikely, we nevertheless wish to carry out a comparison with 
the case for two neutron stars, and use this as a starting point in our investigation. For 
each mass ratio we present the initial conhgurations that were constructed for tidally locked 
binaries, followed by a description of the dynamical runs. 

To investigate the dynamical evolution of the system for a given initial separation, 
we remove the damping term from the equations of motion and give every SPH particle 
and the black hole the azimuthal velocity corresponding to the equilibrium value of hi in 
an inertial frame, with the origin at the center of mass of the system. Each SPH particle 
is assigned a specihc thermal energy Ui=K/{V — 1), and the equation of state is 
changed to that of an ideal gas where P=(P — l)pu. The specihc thermal energy of each 
SPH particle is then evolved individually throughout the simulation according to the hrst 
law of thermodynamics, taking into account the contribution from the viscous terms (see 
Appendix ^). 


5.2. Mass ratio q=l 
5.2.1. Equilibrium configurations 

As described above (section |4.3|) , equilibrium conhgurations were constructed for 
diherent values of the binary separation r. We show in Figure |I| a plot of the total angular 
momentum of the system, J, versus binary separation r. The solid line is the result for two 
point masses in Keplerian orbit. The dotted line is the result of approximating the neutron 
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Fig. 1.— Total angular momentum vs. binary separation for a black hole-neutron star 
binary with mass ratio ^=1. The solid line is the result for two point masses, the dotted 
line is computed using the analytical approach of |Lai, Rasio fc Shapiro 1993b| treating the 
neutron star as a compressible tri-axial ellipsoid, and the crosses are the result of SPH 
relaxation calculations. The Roche limit is at r=2.78. 
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star as a compressible tri-axial ellipsoid ([Lai, Rasio fc Shapiro 1993b|) , and one cross is 
plotted for each SPH calculation at a fixed separation (with iV=8121 SPH particles for 
the neutron star). The presence of a minimum in J is a direct result of the assumption 
of tidal locking of the extended star. Essentially, as the separation is decreased, the spin 
component of angular momentum must grow (just as the orbital component is decreasing), 
and eventually become comparable to the orbital component. The result is that the total 
angular momentum J reaches a minimum and as the separation decreases further, it will 
increase. Such a turning point in the curve of total angular momentum as a function of 
separation marks the onset of an instability, which will lead to orbital decay (|Lai, Rasio & 
^hapiro 1993a|) . Note the strong departure from point-mass behavior, and the excellent 
agreement of the full SPH calculations with the compressible tri-axial ellipsoid treatment 
before the minimum in J is achieved. Close to the minimum, the ellipsoidal approximation 
breaks down and a full numerical treatment is necessary. The cross with the smallest 
value of r corresponds to Roche-Lobe overflow (r=riji=2.78). Our initial determination 
of the stability limit ( Lee fc Kluzniak 1995 ) was higher than this value because our initial 
conditions did not correspond to proper tidal locking. We show in Figure || plots of particle 
positions projected onto the equatorial plane for different separations. As the separation is 
decreased, the effect of the tidal potential becomes more apparent and the neutron star is 
distorted. The conhguration at r=3.0 can be approximated by a tri-axial ellipsoid, while 
for those at r=2.8 and r=2.78 this is increasingly no longer the case. 


5.2.2. Dynamical runs 

We have used several of the conhgurations thus constructed to perform dynamical 
runs. In Figure ^ we show the binary separation as a function of time during several of 
these calculations. It is clear that the conhgurations with initial separations r=3.0 and 
r=2.9 are stable. The separation exhibits oscillations of numerical origin on a period close 
to the orbital period (respectively P=22.81 and P=21.6 in our units, see eq. |^). However, 
the conhguration with r=2.78 and initial period P=20.09 (corresponding to 37 km and 
2.39 ms for P=13.4 km and M=1.4Mq) is unstable, and the separation decays on an orbital 
timescale. We hnd the dynamical stability limit to be rdyn=‘2.78=rjiL. The two diherent 
lines with initial separation r=2.78 correspond to two diherent resolutions (the solid line 
used A^=8121 particles and the dotted line used A^=16944 particles). In the following, we 
will describe the results for our higher resolution run, with A^=16944. In Figure we 
show density contours at various times during the simulation, which was run from t=0 to 
t=100. The left column shows contours in the orbital plane, and the right column shows 
contours in the meridional plane containing the black hole and the center of mass of the 
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Fig. 3.— Binary separation as a function of time for dynamical calculations starting at three 
different separations (r=3.0, r=2.9 and r=2.78) for the black hole-neutron star binary with 
mass ratio q=l. For r=2.78, the solid line used iV=8121 particles and the dotted line used 
A^= 16944 particles. 
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Fig. 4.— Density contours at various times during the dynamical coalescence of the black 
hole-neutron star binary with mass ratio q=l and initial separation r=2.78. The logarithmic 
contours are evenly spaced every 0.25 decades, with the lowest one at logp = —3.25. The 
rotation is counterclockwise and the initial orbital period is P=20.09. The left column shows 
density contours in the orbital plane, the right column shows contours in the meridional plane 
containing the binary axis. The black disk represents the black hole. 
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SPH particles. The black disk with radius represents the black hole. Mass transfer from 
the neutron star to the black hole begins almost immediately (within one orbit) through the 
inner Lagrangian point. The accretion stream gradually becomes thicker and begins to wrap 
around the black hole. By f=50 the star has become considerably elongated and a toroidal 
structure has begun to form around the black hole. The accretion stream then essentially 
breaks, mass transfer stops and a stellar core remains in orbit around the black hole (t=60 
through f=70). The result of this coalescence then appears to be a a stable binary with a 
greatly altered mass ratio and separation {qfinai=0.l9 and r~4.5, corresponding to 60 km), 
consisting of a remnant core with mass Mcore=0.307 (corresponding to O. 43 M 0 ) in orbit 
around a black hole with Mbh=^-Q07 (or 2.25Mq). The mass transfer event is very brief, 
as can be seen in Figure ^ where we plot the mass of the black hole and the mass accretion 
rate onto the black hole as a function of time. The peak accretion rate is dMBH/dt=0.0275 
(corresponding to 0.3MQ/ms). The black hole is still accreting matter towards the end of 
our simulation, but the decreased resolution (recall that accretion entails a loss of particles 
in the simulation) does not allow us to determine the distribution of matter accurately 
beyond f=100 (~11 ms). What can be determined is that the toroidal structure that 
forms does not extend to form a halo that engulfs the black hole completely. The region 
directly above the black hole remains devoid of matter down to the limit of our resolution 
of 7.1 X 10“® (or 10~'^Mq) within ~ 5°. The accretion structure around the black hole was 
not previously observed (|Lee fc Kluzniak 1995| ) because of the low resolution (A^=2176 
particles) of our initial runs. 

The nature of the encounter is reflected in the gravitational radiation waveforms, 
presented in Figure ^ (top). The amplitude of the emitted waves initially rises from 
{tqR/M‘^) h=1.55 to (roi?/M^)hmaa;=l-75, reflecting the decrease in separation. The 
subsequent drop in amplitude and decrease in frequency are a result of the short episode 
of mass transfer, and the hnal waveforms reflect the fact that the binary has survived the 
encounter (the new binary orbital period is P ~ 39 in our units, or 4.5 ms). Thus the 
presence of a remnant core in orbit around the black hole would be immediately apparent 
from an observation of the waveforms emitted during such a coalescence. We also show 
the gravitational wave luminosity in Figure ^ (bottom). For comparison purposes, the 
luminosity for the coalescence of two identical neutron stars is plotted as well, from the 
results of section |4.1| . The maximum luminosity and total energy radiated in gravitational 
waves are {R/MfL^l-^^=0.?,0 (or 1.17 x lO^® erg s-^), ^ 10 (or 

3.4 X 10^2 grg) and {R/MfLg^-^^=0.15 (or 4.5 x 10^^ erg s"^), ^ 6 

(or 2 X 10^^ erg) for the double neutron star and black hole-neutron star case respectively. 
Although in both cases there is still a clear rise and decay in the luminosity curve, the peak 
is much broader ~ 25) than for two neutron stars ~ 15). The slower 






Fig. 5.— Black hole mass (dotted line) and mass accretion rate onto the black hole (solid 
line) as a function of time for the black hole-neutron star binary with mass ratio q=l and 
initial separation r=2.78. 
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Fig. 6.— Top; Gravitational radiation waveforms for the coalescence of a black hole- 
neutron star binary with mass ratio q=l. Two polarizations (h+, left panel; hx, right 
panel) in geometrized units (such that G=c=l) are shown. Bottom: Gravitational radiation 
luminosity in geometrized units (G=c=l) for the coalescence of a black hole-neutron star 
binary (solid line) and a double neutron star binary (dotted line) with mass ratio q=l. The 
constant Lo=c^/G=3.6 x 10^® erg s“^. 
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rise in luminosity and longer timescale of the peak in the black hole-neutron star case are 
a direct consequence of having only one star being deformed and disrupted through tidal 
interactions. 

Finally, we show in Figure the binary separation as a function of time for the 
simulation we have just presented, together with a plot of the separation for a point-mass 
binary with the same mass ratio and initial separation decaying via the emission of 
gravitational waves in the quadrupole approximation. This last curve is given by 
(e.g. [Shapiro fc Teukolsky 1983|) r = r* (1 — , where to = {5rfc^)/{256fiM^G^) and 

Ti is the initial separation at t=0. The total mass is Mt=MBH+M and the reduced mass 
is ij,=MbhM/M f. Clearly, hydrodynamical effects and the instability they induce drive the 
coalescence process on a time scale that is shorter than that purely due to gravitational 
radiation losses, and so we may safely neglect the effect of radiation reaction on the 
evolution of the system over the time period which we have modeled. 

In the calculation described above, we have integrated the thermal energy equation 
according to the hrst law of thermodynamics by setting 

AC = AIT + AQ (4) 

where AQ is calculated using the contribution from the artihcal viscosity (see equation 
in Appendix ^). We will denote this as Case I. To investigate the effect of radiative losses 
on the outcome of this conhguration, we have performed an additional calculation for the 
extreme case in which all the energy produced by viscosity is lost by the system (e.g. to 
neutrinos). This will be denoted as Case II, and is achieved by simply removing the term 
associated with viscous heating from equation (^), so that we now have Af/=AIT. Density 
contours for Case II (with A=16944 particles) are shown in Figure ^ for f=60 and f=80. 
An accretion torus is again clearly visible, and comparing with the contours in Figure ^ for 
Case I, the ring is clearly more conhned both in the orbital plane and along the z-axis. This 
is clearly a consequence of a lack of pressure support due to the removal of thermal energy 
from the system. In Figure ^ we show the rate of energy loss, due to viscous heating. 
The peak rate is {R/M)^Ly ^ 0.65, corresponding to ~ 2 x 10^® erg s“^ and the total 
energy loss from f=0 to t=100 is ~ 5 x 10^^ erg. The mass transfer rate and gravitational 
radiation waveforms are essentially the same as for Case I. There is once more a baryon-free 
axis along the rotation axis of the binary, only this time it is clear of matter within ~ 20 ° 
(down to 4 X 10“®, equivalent to 6 x 16~^Mq) since the torus is more nearly conhned to the 
orbital plane. 


5.3. Mass ratio g=0.8 





Fig. 7.— Binary separation for the coalescence of a black hole-neutron star binary with 
mass ratio q=l (solid line) and for a point-mass binary with the same mass ratio and initial 
separation decaying via the emission of gravitational waves in the quadrupole approximation 
(dotted line). 
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Fig. 8.— Density contours for the dynamical coalescence of the black hole-neutron star 
binary with mass ratio q=l and initial separation r=2.78, Case II. There are eleven evenly 
spaced logarithmic contours down from logp = —0.75 every 0.25 decades. The rotation is 
counterclockwise and the initial orbital period is P=20.09. The left column shows density 
contours in the orbital plane, the right column shows contours in the meridional plane 
containing the binary axis. The black disk represents the black hole. 







. 9.— The thermal energy loss rate due to viscosity in geometrized unit 
c=l'). normalized to Ln=c^/G=3.6 x 10^® erg s“^, as a function of time f 

nd initial seoaration r 
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5.3.1. Equilibrium configurations 


For a mass ratio g=0.8, the black hole mass is M^j:^=1.25 (or 1.75Mq). Hydrodynamical 
effects are still important in the construction of equilibrium configurations in this case, 
as can be seen in Figure where we have plotted the equilibrium angular momentum 
values for a range of binary separations. Approximating the neutron star as a compressible 
tri-axial ellipsoid gives results which are close to the full numerical calculations before the 
minimum in angular momentum is reached, but for r < 3.0 this is no longer the case. We 
identify the separation with minimum angular momentum at r=2.98, and Roche-Lobe 
overflow occurs at rpiL=2.935. For this value of the mass ratio A^=8121 particles model the 
neutron star. 


5.3.2. Dynamical runs 

As for the case with q=l, several equilibrium configurations were used as initial 
conditions for dynamical runs. Figure 0 shows the binary separation as a function of time 
for these calculations. The initial configuration with separation r=2.935 becomes unstable 
on a dynamical time scale. However, it is apparent that the development of the instability 
(and hence the orbital decay) proceeds at a slower rate than for a higher mass ratio. The 
reason for this is that since the black hole is more massive, Roche-Lobe overflow occurs at 
a greater separation than previously observed (for g=l), and so the tidal effects are less 
pronounced on the neutron star. Density contours in the orbital plane for the binary with 
initial separation r=2.935 are shown in Figure Mass transfer proceeds through the 
formation of an accretion stream from the neutron star to the black hole, lasting from f~45 
to f~75 (about 3.5 ms), with a peak accretion rate of 0.017 (0.2MQ/ms). A substantial 
amount of mass is stripped from the neutron star during the encounter (see Figure |^, and 
the final black hole mass is Mbh=^-Q5 (equivalent to 2.3Mq). The core of the neutron star, 
with a final mass Mcore=0.595 (corresponding to O. 83 M 0 ) again survives the encounter. 
However, the most striking difference between this case and the one presented previously 
(g=l) is the absence of a toroidal accretion structure around the black hole. All matter 
stripped from the neutron star is directly accreted by the black hole, and thus in this 
scenario as well the rotation axis above the black hole is devoid of matter. As a result of 
mass transfer the binary separation increases to r~3.6 and thus the final configuration in 
this case is again a stable binary. 

As before, the hydrodynamical evolution of the binary is reflected in the gravitational 
radiation waveforms, presented in Figure |l^ (top). There is an initial rise in amplitude in 
the waveforms, followed by a rapid decline and increase in the period, reflecting the episode 
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Fig. 10.— Total angular momentum vs. binary separation for a black hole-neutron star 
binary with mass ratio ^=0.8. The solid line is the result of the analytical approach of |Lai, 


[Rasio fc Shapiro 1993b| treating the neutron star as a compressible tri-axial ellipsoid, and 
the crosses are the result of SPH relaxation calculations. The Roche limit is at r=2.935. 
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Fig. 12.— Density contours in the orbital plane for a black hole-neutron star coalescence 
with a mass ratio g=0.8 and initial separation r=2.935. The logarithmic contours are evenly 
spaced every 0.25 decades, with the lowest one at logp = —2.75. The black disk represents 
the black hole. Orbital rotation is counterclockwise and the initial orbital period is P=20.58. 
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Fig. 14.— Top: Gravitational radiation waveforms for the coalescence of a black hole- 
neutron star binary with mass ratio g=0.8. Two polarizations (h_|_, left panel; hx, right 
panel) in geometrized units (such that G=c=l) are shown. Bottom: Gravitational radiation 
luminosity in geometrized units (G=c=l) for the coalescence of a black hole-neutron star 
binary with mass ratio g=0.8. The constant Lo=c®/G=3.6 x 10^® erg s“^. 
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of mass transfer and snbseqnent increase in binary separation. The peak amplitnde in this 
case is (roi?/M^)/imaa:=l-93, and again since the binary snrvives the enconnter, the hnal 
amplitnde is not zero bnt (roi?/M^)hjj„aj=1.25. The Inminosity emitted in gravitational 
waves is shown in Fignre 0 (bottom). The general shape of the cnrve resembles very mnch 
that for g=l, as one can expect after stndying the density contonr plots for each case 
(since the toroidal strnctnre that was formed for q=l is essentially azimnthally symmetric 
it does not contribnte signihcantly to the emission of gravitational waves). The cnrve 
presents again a broad peak, with a dnration At~30 (3.4 ms). The peak Inminosity is 
(/?/M)®Lmax=0.16 (or 4.8 x 10^^ erg s“^) and the total energy radiated away in gravitational 
waves is (or 3.4 x 10^^ erg). These valnes are higher than for the 

case with q=l becanse we have normalized onr resnlts to the mass of nentron star, and 
lowering the mass ratio implies the system has a greater total mass. 

In Fignre we compare (as for q=l in Fignre ^ the orbital decay of this binary with 
that for a point-mass binary emitting gravitational waves in the qnadrnpole approximation. 
It is apparent that hydrodynamical effects are less important than for q=l in determining 
the initial orbital evolntion, bnt are nevertheless snbstantial. 


5.4. Mass ratio g=0.31 


Continning the trend towards lower mass ratios, we chose a valne of g=0.31, giving 
^bh=3.22 (corresponding to 4.51Mq) to stndy the evolntion of a binary system where 
most of the mass is contained in the black hole. We again nsed iV=8121 particles to model 
the nentron star. A qnalitative difference appears in the plot of total angnlar momentnm 
as a fnnction of binary separation (Fignre |^. At a separation of r=3.76 (or 50.5 km), 
Roche-Lobe overflow occnrs, so that we cannot constrnct a sequence with decreasing 
separation and constant mass ratio anymore, but this point is reached before the minimum 
in total angular momentum is attained, contrary to what occurred for q=l and g=0.8. Thus 
we expect all conhgurations down to the Roche limit riji=3.76 to be dynamically stable 
(see the discussion for q=l in Section |5.2.1|) . 


The binary separation for a dynamical run with an initial separation r=rpiL=3.76 is 
shown in Figure p!7|. The initial orbital period is P=21.92 (or 2.51 ms). There is no longer 
orbital decay, but rather a slow increase in separation. This is simply due to the fact that 
there is slow Roche-Lobe overflow from the neutron star onto the black hole, transferring 
AM=0.0104 (or O.OldbM©) over ~4.5 orbital periods. The neutron star is not disrupted 
at all during this process, as can be seen in Figure where we plot density contours 
in the orbital plane over the course of the simulation, and in Figure O, where we plot 
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Fig. 15.— Binary separation for the coalescence of a black hole-neutron star binary with 
mass ratio g=0.8 (solid line) and for a point-mass binary with the same mass ratio and initial 
separation decaying via the emission of gravitational waves in the quadrupole approximation 
(dotted line). 
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Fig. 16.— Total angular momentum vs. binary separation for a black hole-neutron star 
binary with mass ratio g=0.31. The solid line is the result of the analytical approach ofE al] 


[Rasio fc Shapiro 1993b| treating the neutron star as a compressible tri-axial ellipsoid, and 
the crosses are the result of SPH relaxation calculations. The Roche limit is at r=3.76. Note 
that Roche-Lobe overflow is reached before the minimum in angular momentum. 
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a) t = 0 b) t=40 



Fig. 18.— Density contours at various times during the dynamical coalescence of the 
black hole-neutron star binary with mass ratio g=0.31 and initial separation r=3.76. The 
logarithmic contours are evenly spaced every 0.25 decades, with the lowest one at log p=-3.0. 
The rotation is counterclockwise and the initial orbital period is P=21.92. The black disk 
represents the black hole. 
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the gravitational radiation waveform for this case. In fact, mass transfer is essentially 
conservative (i.e. the total angular momentum is conserved), and this approximation alone 
can accurately account for the change in separation. To see this, consider the following. 
Conservative mass transfer in a point-mass binary satishes J = = constant, where a 

is the separation and /r is the reduced mass. Since the total mass is conserved and oc 
for a Keplerian binary, we hnd that = constant. For the black hole-neutron 

star binary with initial mass ratio g=0.31, we hnd that 


(^MB Hinitial 

(MbhM) final 


1.01453 


and 


= 1.0167. 

(^initial 

This is consistent with having conservative mass transfer despite the fact that the black 
hole-neutron star binary is not Keplerian. Angular momentum is practically constant since 
there is very little accretion, and this represents the only loss of angular momentum (via a 
spinning up of the black hole, which we do not model). 

This last case is a clear indication of the limitations of the present method of studying 
the coalescence of a black hole with a neutron star. In the hrst place, since the orbit does 
not decay due to hydrodynamical effects, gravitational radiation reaction (which would 
drive the binary mentioned above to coalesce in a time to=31.35, correponding to 3.59 ms) 
cannot be ignored. Second, as is apparent in the density contour plots shown in Figure [^, 
the black hole has become almost as large as the neutron star itself (the Schwarschild radius 
is rs=0.986, or 13.2 km, at the start of the simulation), making the radius of the innermost 
circular stable orbit (ISCO) r/5C'o=3rs=2.958 (or 39.75 km). This means a substantial 
fraction of the neutron star is already at a distance r < risco from the black hole. Thus 
including relativistic effects in this case becomes necessary in order to obtain meaningful 
results. 


6. Discussion 

We present in Table ^the binary separations tjil and corresponding to the Roche 
limit and to the equilibrium conhguration with minimum total angular momentum J for 
each mass ratio. The lack of an entry for when g=0.31 reflects the fact that the orbits 
are stable when the separation is greater than that required for Roche-Lobe overflow. 
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In Table ^ we present a snmniary for the three conhgurations we explored in this 
paper, as well as for the coalescence of two identical nentron stars. The hrst colnmn 
indicates if the coalescence involves two nentron stars or one neutron star and one black 
hole; the second and third columns show the initial and hnal mass ratios respectively; the 
fourth column indicates whether an accretion torus was formed; the hfth, sixth and seventh 
columns show the maximum and hnal gravitational radiation amplitude and the maximum 
gravitational radiation luminosity respectively, all in geometrized units such that G=c=l; 
the eighth column shows the mass of the central object for the case of the double neutron 
star coalescence, and that of the remnant core in the case of the black hole-neutron star 
coalescence. 


Note that only in the simulation with a high mass ratio ((3'=1) did an accretion 
structure appear around the black hole. Otherwise, all matter stripped from the neutron 
star was directly accreted by the black hole. A striking result is that in every case, the 
neutron star is not completely disrupted, but a remnant core survives in orbit around the 
more massive black hole. Due to the violent mass transfer episode (for ^=0.8 and q=l), this 
core is transferred to a higher orbit on a time scale of one orbit. For now, we are unable 
to model the system on longer time scales than we have presented here, but presumably 
the emission of gravitational waves will make the orbit decay by decreasing the separation. 
Clearly mass transfer will occur before the separation reaches zero, as soon as the neutron 
star overflows its Roche Lobe. For the equation of state we have considered (with F=3), 
the mass-radius relationship is i? oc and one can estimate the size of the Roche Lobe 
as Rj?L=0.46a(M/Mi)^/^ ([Paczyhski 1967|) . Using the rate of angular momenutum loss 
calculated from the quadrupole approximation, given by (phapiro fc Teukolsky 198^) 


dJ _ 32 

dt 5 ’ 

we can hnd the separation at which mass transfer will begin anew, as well as the time it will 
take for the orbit to shrink to this size. Applying this reasoning to the binary with initial 
q=l after the initial mass transfer event, we hnd that a second episode of mass transfer will 
occur ~40 ms later, and for the case with initial g=0.8, the delay will be ~4 ms. Thus the 
length of the coalescence process is extended from a few milliseconds to possibly several 
tens of milliseconds. 


Another important result is the fact that in every case, there is a baryon-free line of 
sight to the black hole along the rotation axis of the binary throughout the simulation. As 
stated in section the limit we infer for the amount of matter along this axis from 

our calculations is within about 5° of the rotation axis. It is set by the resolution 

in our calculation and thus represents an upper bound at this point. For the calculations 
with a lower mass ratio (g=0.8 and g=0.31) there was no visible accretion structure present 











38 


around the black hole, and so the rotation axis was clear of matter as well. These last two 
results may make the coalescence of a black hole with a neutron star a promising candidate 
source for the production of GRBs ([Kluzniak &: Lee 1998|). 


The survival of the neutron star core suggests a different outcome than what was 
initially expected for the coalescence of a black hole with a neutron star ([Wheeler 1971|) , 
and identihes this process as the only one known capable of producing low-mass neutron 
stars. The violent episode of mass transfer strips the neutron star of so much matter that it 
may drive it below the minimum mass required for stability. If this happens, then the core 
explodes ( [Blinnikov et al. 1984| , [Sumiyoshi et al. 1998[ , [Colpi, Shapiro fc leukolsky 1991[ ) in 
approximately 0.1 seconds after a slow expansion phase that may last 20 seconds ([Sumiyoshi 


\et al. 19^5[). The presence of the black hole would complicate the outcome of this event. 


This work was supported in part through KBN grant 2P03D01311 and DGDPA-UNAM. 


A. Numerical Code 

The equations of motion in SPH are essentially those of an N-body problem, with each 
particle representing a fluid element. Thus the equations for particle i can be written as: 

hi = Vi, 


rriiVi — Fic+FiH^ 

where r^, Vi and m* are the position, velocity and mass of particle i respectively, and FiG 
and FiH denote the gravitational and hydrodynamic forces. The density at the position of 
particle i is 


Pi = p{.ri) = ^rrijWij, 

3 

where Wij is the symmetrized kernel for particles i and j. With hi representing the 
smoothing length for particle i, we have 

Ihij ^^(1 ^3 I) hp), hij —(^hi T hj). 

It is convenient to perform calculations with a kernel that has compact support. We 
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use the form of Monaghan fc Lattanzio (1985) 


hh(r, h) = 


Tih^ 


i-lh) 

i(2 


+ 




4 

0 , 


)' 


0 < r/h < 1, 

1 < r/h < 2, 

2 < (r/h). 


To obtain a value of the pressure for any particle, we use the equation of state. For an 
ideal gas this becomes 

Pi = {T - l)uipi, 

where F is the adiabatic index and Ui is the thermal energy per unit mass for particle i. 


FiH includes the contribution from the pressure gradient and from artihcial viscosity. 
In symmetrized form, it can be written as 


where 


P^Pj 


FiH = - ^rrijmi ( 2- -^ + 11 *^ | ViWij, 

PiPj 


^ f {-acijpij + Pip Vij ■ Tij < 0 

1 0 , Vij ■ Vij > 0 ’ 


(Al) 


Pij = 2 , 2 1.2 • 

Pij + h %■ 

Here Vij=Vi—Vj, rij=ri—rj, a and /3 are constants, q is the speed of sound at the position 
of particle i, and Cij={ci + Cj)/2 and Pij={pi + Pj)/2. In the simulations presented in this 
paper we have used a=l, (3=2 and 


The form of FiH , and in particular the symmetrization of the pressure gradient term, 
ensures that forces between pairs of particles are symmetric, and so linear and angular 
momentum are automatically conserved. The viscosity is introduced with a term linear in 
the velocity differences, which produces shear and bulk viscosity, and a quadratic term to 
handle high Mach numbers in shocks, analogous to the Von Neumann-Richtmyer artihcial 
viscosity used in hnite-difference methods ([Richtmyer 1967|) . 


The value Ui of the thermal energy per unit mass is evolved in time according to the 
the hrst law of thermodynamics, taking into account the contribution from the viscous 
terms. The symmetrized form is: 


dui 

dt 



+ Ifj 


Vij • ^7iWij. 


(A2) 
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The thermal energy equation is integrated by essentially using a two-step 
procedure (|Hernquist fc Katz 1989| ) in order to preserve accuracy. This is necessary 
because the heating rate depends on the pressure, and via the equation of state, on the 
thermal energy itself. 


Individual smoothing lengths hi are adjusted so that every particle has an approximately 
constant number v of hydrodynamical neighbors, i.e., those for which rijjhij < 2. This 
ensures that adequate spatial resolution is maintained and that the level of accuracy of the 
interpolations remains uniform throughout the fluid. If hi^n is the smoothing length for 
particle i at step n, the value at step n + 1 is found according to (|Hernquist fc Katz 19891) : 


ft”+‘ 



(A3) 


where z/"" is the number of hydrodynamical neighbors at step n. We have found that a value 
of z/=64 yields adequate sampling of the fluid without requiring excessive computation time. 

Gravitational force calculations constitute the most time-consuming part of an iV-body 
code. A direct calculation requires 0{N^) operations, and this becomes prohibitive as the 
number of particles is increased above a few hundred. In order to make this process more 
efficient we have incorporated a binary tree structure (|Hernquist &: Katz 1989|) into our code 
to obtain gravitational forces. This makes the number of required operations of O(A^IogA^), 
thus allowing for simulations with several thousand particles. The tree structure is built by 
pairing two particles into a new pseudo-particle, called a node, if they are mutual nearest 
neighbors, and proceeding with this pairing over all particles and newly created nodes 
until there is only one node, at the center of mass of the particles being considered. This 
hierarchical association makes up the binary tree which we use to compute the forces (for a 
detailed description see |Benz et al. 199Ci|) . The force on any particle can then be obtained 
by adding up contributions from particles close enough to be considered individually and by 
more distant aggregations, which can be described in terms of their multipolar expansion 
up to the quadrupole term. To decide if an aggregation of particles must be resolved to 
obtain the force on particle i, we consider the ratio s/d, where d is the distance from 
particle i to the aggregate and s is its size. If s/d > 9, where d is a dimensionless tolerance 
parameter, the aggregate is resolved further. If the converse is true, the contribution to 
the force is obtained from the quadrupole expansion. We adopt 6*=0.8, which provides an 
accuracy for the force of better than one part in 10^. However, the aggregate is resolved 
further if it contains any individual particles which may be hydrodynamical neighbors for 
particle i. The subroutine corresponding to this approach simultaneously provides F^g and 
the hydrodynamical neighbors for every particle that are used to obtain FiH and dui/dt 
in equations ([Al|) and (|A2|) . Our calculations of the gravitational forces are completely 
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Newtonian. 

We use a leapfrog algorithm accurate to second order to integrate the equations of 
motion. The time step is adaptive, and is taken to satisfy a combination of the Courant 
criterion for stability and a further restriction on the maximum change allowed in velocity 
for any particle during one time step to conserve accuracy. We follow [Monaghan (1989l| and 
set the time step to be At=min {Ati, At 2 ) with 

Ati = min(hj/hj)^'^^ 


At 2 = 0.15 min 


hi 


^Ci + 1.2acj + 1.2/9 maxj(pjj) 


One of the main objectives of the present work is to investigate the gravitational waves 
emitted by coalescing binaries. In our Newtonian simulations we calculate the waveforms 
of the emitted waves and the total luminosity in the quadrupole approximation. The 
amplitude of the retarded wave at time t and a distance tq away from the source is given 
by (jMisner, Thorne fc Wheeler 1970|) 


h 


TT 

jk 


2 G ■■ TT 

Tq G 



1 


where Ijk Pjn^nmPmk PjkPlmhlm/2: 
is the transverse traceless part of the reduced quadrupole moment, 

hjk hjk 


with Ijk = 



and J = ^ Ijj. 
j 


The Pij are the components of the projection operators onto the plane transverse to 
the radial direction: P^j = Sij — fiiUj with rii = Xi/rQ. 

The luminosity of the gravitational waves at a distance r from the source is 


_ dE 
dt 


IG 
5 c® 


hdjk ddji^ 


(A4) 


dt dt 
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We use the method of [hinn (1989) and obtain the formula applicable to SPH for the 
second derivative of the quadrupole moment ( [Rasio fc Shapiro 1992| ): 


'fjk _ 

fluid 


+ 


Pi 




The summation is over the SPH particles and the superscripts indicate the cartesian 
components. We identify the gravitational acceleration on particle i as gi. The term 
only includes the contribution from the SPH particles contained in the neutron star. We 
add the terms arising from the presence of the black hole as a point mass contribution: 

P^H = Mbh (fpBH'^BH + ^Bh9bH + Pbh9%h) i 


to obtain the total value of I as 


pk 



I fjk 
^ ^fluid- 


The gravitational waveforms are thus obtained directly from dynamical and 
hydrodynamical variables of the system, and one numerical differentiation of is required 
to obtain the luminosity. 


The amplitudes of the waveforms for an observer a distance tq away along the axis of 
rotation of the binary are given by 


roK 



(AS) 


rohx 


G ■■ 

2 _ f 


(A6) 
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Table 1. Comparison of results for the coalescence of two identical polytropes with a stiff 

equation of state (T = 3). 



{roR/ M^)h^ax'" 

{roR/M^)h final"" 

{R/Mf{Ln,ax/Lo)^ 

Pc 

Mhalo/Mtot 

This work 

2.2 

0.2 

0.39 

0.44 

0.16 

RS 

2.2 

0.2 

0.37 

0.40 

0.18 


^The peak and final gravitational radiation amplitudes at a distance tq from the source 
(see equations ||A5|| and [|A6||). 


'’See equation (|X^). Lq = /G = 3.59 x 10^® erg s ^ 


Table 2. Critical Separations for black hole-neutron star binaries. 


Q ^ min 

trl 

1.00 2.82 

2.78 

0.80 2.97 

2.94 

0.31 ■■ ■ 

3.76 
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Table 3. Summary of results for the coalescence of two neutron stars and a black hole 

with a neutron star? 



Qi 

Qf 

Torus 

{roR/ M‘^)hmax 

(roR/M^)h final 

{R/Mf{Ln,ax/Lo) 

^core 

NS-NS 

1.000 

1.000 

no 

2.20 

0.20 

0.390 

0.82 

BH-NS 

1.000 

0.190 

yes 

1.75 

0.55 

0.145 

0.3070 

BH-NS 

0.800 

0.360 

no 

1.93 

1.25 

0.160 

0.5950 

BH-NS 

0.310 

0.306 

no 

3.59 

3.48 

0.425 

0.9896 


^Notation is the same as in Table 






